#
# Copyright (c) 2014 Luciano Xumerle. All rights reserved.
#
# This file is part of gephinfo.
#
# gephinfo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

package trimSEQ;

require 5.003;
use warnings;
use strict;
use utf8;

=head2 new( id, flag, chrom, pos, mapq, cigar, chromMate, posMate, templateLength, seq, qual )

The contructor of the class.
Input is simply a splitted SAM line

 1 QNAME String [!-?A-~]{1,255} Query template NAME
 2 FLAG Int [0,216-1] bitwise FLAG
 3 RNAME String \*|[!-()+-<>-~][!-~]* Reference sequence NAME
 4 POS Int [0,231-1] 1-based leftmost mapping POSition
 5 MAPQ Int [0,28 [ -1] MAPping Quality
 6 CIGAR String \*|([0-9]+[MIDNSHPX=])+ CIGAR string
 7 RNEXT String \*|=|[!-()+-<>-~][!-~]* Ref. name of the mate/next read
 8 PNEXT Int [0,231-1] Position of the mate/next read
 9 TLEN Int [-231+1,231-1] observed Template LENgth
 10 SEQ String \*|[A-Za-z=.]+ segment SEQuence
 11 QUAL String [!-~]+ ASCII of Phred-scaled base QUALity+33

=cut

sub new
{
    my $class          = shift;
    my $id             = shift;
    my $flag           = shift;    # 0 in single mate forward alignment
    my $chrom          = shift;
    my $pos            = shift;
    my $mapq           = shift;
    my $cigar          = shift;
    my $chrommate      = shift;    # * in single mate
    my $posmate        = shift;    # 0 in single mate
    my $templatelength = shift;    # 0 in single mate
    my $seq            = shift;
    my $qual           = shift;

    #
    # INITIALIZE CLASS
    #
    my $self = {
        'ID'             => $id,
        'FLAG'           => $flag,
        'CHROM'          => $chrom,
        'POS'            => $pos,
        'MAPQUALITY'     => $mapq,
        'CIGAR'          => &SAMtrim::cigar2list($cigar),
        'CHROMMATE'      => $chrommate,
        'POSMATE'        => $posmate,
        'TEMPLATELENGTH' => $templatelength,
        'SEQUENCE'       => $seq,
        'QUALITY'        => $qual,
    };
    bless( $self, $class );
    return ($self);
}

=head2 getID()

Return the read ID

=cut

sub getID
{
    my $self = shift;
    return $self->{'ID'};
}

=head2 getChrom()

Return the read chromosome

=cut

sub getChrom
{
    my $self = shift;
    return $self->{'CHROM'};
}

=head2 getMapQuality()

Return the mapping quality

=cut

sub getMapQuality
{
    my $self = shift;
    return $self->{'MAPQUALITY'};
}

=head2 getCigar()

Return the current cigar

=cut

sub getCigar
{
    my $self = shift;
    return join '', @{ $self->{'CIGAR'} };
}

=head2 cigar2list( $cigar )

Static method to split the cigar element.
The method is used by the class constructor.

=cut

sub cigar2list
{
    my $cigar = shift;
    return ['*'] if ( $cigar eq '*' );
    my $res   = [];
    my $val   = '';
    my @cigar = split '', $cigar;
    while (@cigar)
    {
        my $e = shift @cigar;
        if ( $e =~ m/\d/ )
        {
            $val .= $e;
        }
        else
        {
            next if ( $e eq 'H' );
            push @$res, $val;
            push @$res, $e;
            $val = '';
        }
    }
    return $res;
}

=head2 readInfo()

Returns the list of four coordinates:
[ start, stop, softClippedStart, softClippedStop ]

 CIGAR FIELDS
 ^^^^^^^^^^^^
 M 0 alignment match (can be a sequence match or mismatch)
 I 1 insertion to the reference
 D 2 deletion from the reference
 N 3 skipped region from the reference
 S 4 soft clipping (clipped sequences present in SEQ)
 H 5 hard clipping (clipped sequences NOT present in SEQ)
 P 6 padding (silent deletion from padded reference)
 = 7 sequence match
 X 8 sequence mismatch

=cut

sub readInfo
{
    my $self = shift;
    return [ 0, 0, 0, 0 ] if ( $self->{'CIGAR'}->[0] eq '*' );

    my $start = $self->{'POS'};
    my @cigar = @{ $self->{'CIGAR'} };

    my $clipstart = $start;

    $clipstart = $start - $cigar[0] if ( $cigar[1] eq 'S' );
    my $noclipend = $clipstart;
    my $end       = $clipstart;
    for ( my $i = 0 ; $i < $#cigar ; $i += 2 )
    {
        if ( $cigar[ $i + 1 ] ne 'I' )
        {
            $end += $cigar[$i];
        }
    }
    $end--;
    $noclipend = $end;
    $noclipend = $end - $cigar[ $#cigar - 1 ] if ( $cigar[$#cigar] eq 'S' );
    return [ $start, $noclipend, $clipstart, $end ];
}

=head2 trimSoftClip()

Removes the soft clipping from start and end of the reads
and returns the list
[ trimmedSequence, trimmedQuality, [cigar list] ]
The method updates the internal values stored into the class object.

=cut

sub trimSoftClip
{
    my $self = shift;
    return [ $self->{'SEQUENCE'}, $self->{'QUALITY'}, $self->{'CIGAR'} ]
      if ( $self->{'CIGAR'}->[0] eq '*' );
    my $seq   = $self->{'SEQUENCE'};
    my $qual  = $self->{'QUALITY'};
    my @cigar = @{ $self->{'CIGAR'} };

    if ( $cigar[$#cigar] eq 'S' )
    {
        $seq  = substr( $seq,  0, -$cigar[ $#cigar - 1 ] );
        $qual = substr( $qual, 0, -$cigar[ $#cigar - 1 ] );
        pop @cigar;
        pop @cigar;
    }

    if ( $cigar[1] eq 'S' )
    {
        $seq  = substr( $seq,  $cigar[0] );
        $qual = substr( $qual, $cigar[0] );
        shift @cigar;
        shift @cigar;
    }

    $self->{'SEQUENCE'} = $seq;
    $self->{'QUALITY'}  = $qual;
    $self->{'CIGAR'}    = \@cigar;

    return [ $seq, $qual, \@cigar ];
}

=head2 reverseCigar( \@cigarList )

Static method used to reverse a cigar list.
Input is a cigar list.
Output is the pointer to the reversed cigar list.

=cut

sub reverseCigar
{
    my @cigar = @{ (shift) };
    return ['*'] if ( $cigar[0] eq '*' );
    my $reversed = [];

    for ( my $i = 0 ; $i < $#cigar ; $i += 2 )
    {
        unshift @$reversed, $cigar[ $i + 1 ];
        unshift @$reversed, $cigar[$i];
    }
    return $reversed;
}

=head2 reverseSequenceQualityCigar()

Reverses the sequence (and complement the bases),
the quality string and the cigar list.
The method updates the class values and returns the
list with the changed elements:
[ trimmedSequence, trimmedQuality, [cigar list] ]

=cut

sub reverseSequenceQualityCigar
{
    my $self = shift;

    my %rev = (
        'T' => 'A',
        'A' => 'T',
        'G' => 'C',
        'C' => 'G',
        'N' => 'N',
        '.' => '.'
    );

    my @ss = split '', scalar reverse $self->{'SEQUENCE'};
    my $seq = '';
    foreach (@ss)
    {
        $seq .= $rev{$_};
    }

    $self->{'SEQUENCE'} = $seq;
    $self->{'QUALITY'}  = scalar reverse $self->{'QUALITY'};
    $self->{'CIGAR'}    = &SAMtrim::reverseCigar( $self->{'CIGAR'} );

    return [ $seq, $self->{'QUALITY'}, $self->{'CIGAR'} ];
}

=head2 mapOneRead( $bedIntervalOBJ, $index )

Checks if the coordinates of the read are
in according with an amplicon stored into
$bedIntervalOBJ and trim the read if longer
than the matched amplicon.

=cut

sub mapOneRead
{
    my $self  = shift;
    my $bed   = shift;
    my $index = shift;

    return if ( $self->getCigar() eq '*' );
    my $p = $self->readInfo();

    #
    # find the possible amplicon coordinates
    # if NO, then return
    #
    my $astart = 0;
    my $aend   = 0;

    my $stopboh = $bed->getStop( $self->getChrom(), $p->[0] );
    $stopboh = $bed->getStop( $self->getChrom(), $p->[2] ) if ( $stopboh < 0 );
    my $startboh = $bed->getStart( $self->getChrom(), $p->[1] );
    $startboh = $bed->getStart( $self->getChrom(), $p->[3] )
      if ( $startboh < 0 );

    #print join "\t", @$p,  $startboh,  $stopboh , "\n";
    if ( $stopboh > -1 )
    {
        $astart = $p->[0];
        $aend   = $stopboh;
    }
    elsif ( $startboh > -1 )
    {
        $astart = $startboh;
        $aend   = $p->[1];
    }
    else
    {
        return if ( $astart == $aend );
    }

    $self->trimSoftClip();

    #
    # tolgo da inizio o fine una base alla reads
    # una base al cigar
    # una base alla qualità
    # una base a start o end della sequanza
    #
    my @seq  = split '', $self->{'SEQUENCE'};
    my @qual = split '', $self->{'QUALITY'};
    my @cigar = @{ $self->{'CIGAR'} };

    my $sstart = $self->{'POS'};

    #
    # ACCORCIA TESTA READ
    #
    while ( @cigar && $astart > $sstart )
    {

        # con insert la posizione sul reference
        # resta uguale e tolgo le basi
        if ( $cigar[1] eq 'I' )
        {
            for ( my $i = 0 ; $i < $cigar[0] ; $i++ )
            {
                shift @seq;
                shift @qual;
            }
            shift @cigar;
            shift @cigar;
        }

        # con delezione non tolgo basi
        # ma mi sposto sul reference
        elsif ( $cigar[1] eq 'D' )
        {
            $sstart = $sstart + $cigar[0];
            shift @cigar;
            shift @cigar;
        }
        else
        {
            $sstart++;
            shift @seq;
            shift @qual;
            if ( $cigar[0] == 1 )
            {
                shift @cigar;
                shift @cigar;
            }
            else
            {
                $cigar[0]--;
            }
        }
    }

    #
    # ACCORCIA CODA READ
    #
    my $stop = $p->[1];
    while ( @cigar && $stop > $aend )
    {
        if ( $cigar[$#cigar] eq 'I' )
        {
            for ( my $i = 0 ; $i < $cigar[ $#cigar - 1 ] ; $i++ )
            {
                pop @seq;
                pop @qual;
            }
            pop @cigar;
            pop @cigar;
            last if ( $stop - 1 <= $aend );
        }
        elsif ( $cigar[$#cigar] eq 'D' )
        {
            $stop = $stop - $cigar[ $#cigar - 1 ];
            pop @cigar;
            pop @cigar;
        }
        else
        {
            if ( $cigar[ $#cigar - 1 ] == 1 )
            {
                pop @cigar;
                pop @cigar;
            }
            else
            {
                $cigar[ $#cigar - 1 ]--;
            }
            $stop--;
            pop @seq;
            pop @qual;
        }
    }

    print join( "\t",
        $self->getID() . '-' . $index, 0,
        $self->getChrom(),             $sstart,
        $self->getMapQuality(), ( join '', @cigar ),
        '*', 0,
        0, ( join '', @seq ),
        ( join '', @qual ) ),
      "\n"
      if ( $#cigar > 0 );
}

=back

=head1 AUTHOR

Luciano Xumerle (luciano.xumerle@gmail.com)
